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Abstract, Psychophysical experiments show that the perceived colour of an object is 
relatively independent of the spectrum of the incident illumination and depends only on the 
surface reflectance. We demonstrate a possible solution to this underdetermined problem 
by expanding the illumination and surface reflectance in terms of a finite number of basis 
functions. This yields a number of nonlinear equations for each colour patch. We show that 
given a sufficient number of surface patches with the same illumination it is possible to solve 
these equations up to an overall scaling factor. Generalizations to the spatial dependent 
situation are discussed. We define a method for detecting material changes and illustrate a 
way of detecting the colour of a material at its boundaries and propagating it inwards. 
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A Method for Computing Spectral Reflectance 


1. introduction. The colour of an object depends on the spectral distributions of the 
incident light and the reflectance functions of the object. More precisely at each point x 
the receptors (for example, the rods and cone of the eye) measure a quantity given by' 

l^x) = / o /1 (X)7?(®, \)S(x, \)d\ (t.l) 

where o M (X) is the absorption function of the receptors (/;. = 1 to 4 if we consider the rods 
and cones), E(x,\) is the incident illumination and S(x,\) the reflectance of the surface. 
In terms of the notation used by Hurlbert (1985) our function S{x,\) is the product of the 
albedo p{x,\) and the reflectance function R(n,k,s). 

The.se quantities l lt (x) are then evaluated and combined to give the perception of colour. 
For the human eye it is assumed that the outputs of the three cones correspond to the 
perception of colour, The three cones corresponds to the three dimensions of the space of 
perceived colours (von Helmholtz 1866, or see Feynman 1965 for an interesting review). 1 

A number of experiments show that perceived colour of an object is relatively independent 
of the spectral distribution of the incident light . This is the colour constancy effect 
demonstrated most clearly by Land on Mondrians. Since the incident illumination and the 
surface reflectance are confounded in the input (1.1) some assumptions must be made to 
disentangle them and to obtain a colour which depends on the surface reflectance alone. 
This paper suggests a scheme for doing this. Alternative schemes are described in Rubin 
and Richards (1984), Hurlbert and Poggio (1985). 

This paper assumes that both the colour, the incident illumination and the surface reflectance 
can be expressed as a finite sum of basis functions in spectral space. 2 It is then shown 
that given a sufficient number of patches of different surface reflectivity with roughly the 
same incident illumination it is possible to solve for the surface reflectance up to an overall 
scaling factor. In such a theory the colour of a patch would depend on the colour of 
the neighbouring patches, this is consistent with Land’s experiments. Rubin and Richards 
(1984) have considered ways of distinguishing such "material" patches by considering the 
spectral properties of the observed intensities. They discuss "lawful processes" (such as 

‘The rods are not normally assumed to contribute to colour perception, however, we will argue later 
that although their inputs do not directly contribute to perceived colour they can still be used to 
"normalize" colour perception. 

“This corresponds to the standard regularization approach (Poggio and Torre, T984) of solving an 
underconstrained problem by restricting the space of solutions. 
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orientation and shadowing) which change the albedo while still corresponding to the same 
material and propose an "opposite sign slope condition" to detect such change. 

The first part of this paper, although done independently, is similar to work of VVandell and 
Maloney (1984a, 1984b) :! . In their work the illuminant and spectral surface reflectances 
are also expanded in terms of a finite number of basis functions in spectral space and 
neighbouring colour patches, with the same assumed illumination. Our formula for the 
number of patches required differs from theirs, however, and furthermore we demonstrate 
a simple way to solve the equations. Our scheme can also be extended using results from 
Hurlbert and Poggio(1985) to cases when the spectral components of the illuminant and the 
reflectance function vary with position. For this case we can find the colour of an object at 
its boundaries and then fill in the interior. This method will automatically find the material 
boundaries and thus will act as a material edge detector in the sense of Rubin and Richards 
(1984). 

2. The Basis Functions. First we consider the case when there is no spatial dependence 
of the illuminant or the surface reflection. The inputs to the receptors are then 

l tl = J aii {\)E{\)S{\)d\. (2.1) 

We now expand E(\) and £(X) in terms of basis functions Z?,(X) and Cj{\) 


n m 

m = £ TiBifr), S(\) = £ 7j cy(x). (2.2) 

‘=1 i—i 

For the present we will set n = m = 3 and choose the same basis functions 

3 3 

E{\) = £ r t -B f (X), S(X) = £ 7^.(X). (2.3) 

i'=i j—i 

This is merely for convenience and because of psychophysical evidence that perceived 
colours form a three dimensional space (Helmholtz, (1866), for a review see; Feynman 
1965). It is straightforward to generalize the analysis to different values of n and m and for 
different numbers of rods and cones, Substituting (2.3) into (2.1) gives 

l M - r nj j <v(X)i3,(X)/ly(X)eiX (2.4) 

where we use the summation convention over i and j. 

W e can w rite this as 

3 I am grateful to Ted Adelson for bringing their work to my attention. 
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l it TfjTi'ij. (2.5) 

The (T£) depend only on the absorption coefficients and the basis functions. They are 
therefore parameters of the system and hence are known (the system might be taught to 
learn the optimal (T'T)). The (r ; ) and (7 3 ) describe the iliuminant and the surface reflection 
respectively. We therefore have four equations (n = 1 to 4) for six unknowns {i = 1 to 3 
, j — 1 to 3). Note, however, that there are effectively only five unknowns because of the 
scaling ambiguity 


Ti i-> ^ (2.6) 

This ambiguity 4 can be traced back to (2.1) where it corresponds to the scaling 

E{\) P E{\), S{\) i ► (2.7) 

This scaling ambiguity cannot be dealt with by our theory and a further normalization 
is required. Note this is the only ambiguity in our theory and only one normalization is 
required. This is in contrast to the theories of Land(1983), Horn(1975) and Blake(1984) 
where colour is computed independently in seperate channels each of which has to be 
normalized seperately. Unlike these theories (for a review of these see Hurlbert 1985) our 
method allows interactions between channels and hence we only need one normalization. 

To resolve this ambiguity we can normalize the light source by making (r ; ) a unit vector, or 
by setting n = 1. 

A single patch gives us four equations for five unknowns. Suppose we have a neighbouring 
patch with a different spectral surface reflectance. This second surface reflectance can be 
expressed in terms of basis functions by 

3 

S*(X) - £ (2.8) 

1=1 

Now if we make the reasonable assumption that the spectral distribution of the iliuminant is 
unchanged on the second patch we measure four new quantities 


■’whose nature was clarified during discussions with Hurlbert and Poggio 


(2.9) 
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which involve only three new unknowns (fy). We combine this with (2,8) to get eight 
equations for eight unknowns. Number counting suggests this is enough to get a unique 
solution. The equations are non-linear however and, as we will show in the next section, in 
some cases there may be more than one solution. 

For the human eye to use a colour scheme like this with four receptors, the rods are 
necessary as well as the cones. If we just use the cones and have three receptors it follows 
from (2.10) that we can only have two basis functions and hence only a two dimensional 
space of perceived colour. Observe that the output from the rods can be treated differently 
from the output from the cones. It is combined with them to solve (2.4) and (2.6) and 
determine the (fy) and the (qy) but it is only the inputs from the cones that gives the 
colours. The rods and the cones combine to determine the perception of colour from the 
spectral information at the cones. 

We now consider the general case when there are p receptors, n basis functions for the 
illuminant and m for the reflectance function. Suppose we have q patches. Then the number 
of equations is pq. The number of variables is n + mq - 1 (the 1 comes from the scaling 
ambiguity). To obtain a unique so!ution 5 we need 


pq > n + mq - 1. (2-10) 

This is a neccessary condition. To check it is sufficient we will show that the non-linear 
equations can be solved. 

Note from (2.10) that we always need p > m, that is the number of different photoreceptors 
must be greater than the number of basis functions for the reflectance. Thus to obtain a 
three dimensional colour space with this scheme at least four photoreceptors are required.) 

The number of patches required is 


q > 


n — 1 
p — m 


( 2 . 11 ) 


So if there are two basis vectors for the illuminant and surface reflectances (n-m = 2) and 
three receptor fields (p = 3) we only need one patch. 

3. Solving the Equations. 

We now consider the equations in more detail. They are of the form 


yT(‘. 


(3.1) 


5 This equation has also been derived by J. Rubin (pers. comm.). It does not imply that g > n, a 
formula stated by Wandell and Maloney (1984a), and a counterexample is given in the next section. 
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The matrices 7’T are parameters of the system given by 

ni = / (3.2) 

If the basis functions (A(X)) cJo not overlap then the off-diagonal terms of T(‘- are zero 
(rfr — 0 , i 7 ^ j). This case is considered by Hurlbert and Poggio (1985). This simplifies the 
equations (3.1) but at the cost of preventing interactions between different colour channels 
(although interactions can be introduced later). These interactions are important for our 
purposes since they mean that only one colour normalization (the overall scaling given by 
(2,7)) is needed. 

The basic strategy to solve the equations (3.1) is to take linear combinations of the matrices 
T'-j to obtain simpler equations. We illustrate this by considering the case when there are 
three photoreceptors, and two basis functions (n — m — 2, p — 3). As shown in the previous 
section only one patch is needed to solve the equations. These are 

l M = miT%, i = 1,2 .j - 1, 2. m = 1,2,3. (3.3) 

The ( 7 ’T) are three two-by-two symmetric matrices. The space of such matrices is three 
dimensional and so the 7'T form a basis for this space . 6 Hence we can find three vectors 
p“ (a — J, 2 ,3) such that 


1 O' 
.0 0 . 



'O' O' 
.0 1 . 



o r 

.1 o. 


E 


p 3 7’f*.. 
* 5 * i] 1 


Combining (3.4) with (3.3) defines three (known) numbers A u Az,A 3 by 


(3.4) 


3 3 3 

A i = P/V/i = nil, A = P 1 1 I> — r '^ 2 > A * — ]C p l l > 1 = Tl T2 + r 2 71 . (3.5) 

n — 1 (i — 1 n—l 

We normalize the light source by setting rj — 1 . This reduces (3.5) to 


__ ’ 7l — A l>l2 T 2 — A ‘i,AiT2 + 72 — A 3 . (3.8) 

6 We choose the (72)) so that they are linearly independent. We are free to do this since they are the 
parameters of the system. 
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These can be solved to give 


. Ai ± (A\ — T (/lf ( — AAiA^Y^ 

n = ) f 2 -- 2 / l ~- ~~~’ 71 = 1,72 ™- 2 - 



This gives two possible solutions (the values of r 2 and 72 are interdependent). The constraint 
that the illuminant and the reflectance function are both positive functions may not help to 
resolve this, the two roots of both r 2 and 72 all have the same sign. 


Mote that only one patch is needed to solve for this case although there are two basis 
functions. 

The General Case 


The basic strategy to solve the equations consists of taking linear combinations of the (l v / 3 ) 
and reducing them to some algebraic equations having only a finite number of solutions. 


For four photoreceptors and three basis functions we require two patches with the same 
illumination. We have two sets of equations, one from each patch, 


/"""N 




Ifj. = r rij n 3 , k, = (3.8) 

A simple way to solve these equations is by choosing the basis vectors (/T(X)) so that the 
(T'l) generate a space spanned by 


"1 

0 

(f 


'0 

0 

O ' 


'0 

1 

O ' 


'0 

0 

o ' 

0 

0 

0 

> 

0 

1 

0 

> 

1 

0 

t 

) 

0 

0 

0 

.0 

0 

0 - 


.0 

0 

0 . 


.0 

1 

0 . 


.0 

0 

1 . 


The conditions on the basis vectors (/? t (X)) which this requires are 

J a M j5i(X)£ta(X)<fX — a' 1 Bi{\)B 2 (\)d\ — J a> 1 B 2 {\)B 3 {\)d\. (3.10) 

To see this observe that the space of symmetric three by three matrices is six dimensional. 
The matrices defined by (3.9) define a four dimensional subspace and to ensure that the 
(rfj) lie in this subspace we must impose the two conditions of (3.10). To check the 
exact form of (3.10) consider all possible matrices i7 ti that can be generated by linear 
combinations of (3.9). First observe that it is impossible to get a matrix with II 13 =/ 0, this 
gives us the first equation of (3.10). The second equation comes from noting that we must 
have II 12 — H 32 for all II. 

We now take linear combinations of (3.8) using four vectors 0 ® (a — 1,2,3,4) as in (3.4) and 
(3.5) until we generate the basis matrices of (3.9). The equations (3.8) are now of form 
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(3.11) 

(3.12) 
and 

A 1 (t 2 2 t,) + 7l 2 (r 3 + t\) + A 3 r! = A 4 t 2 t 3> B 1 (t 2 t 3 ) + 0 2 (r 3 + r 2 ) + B 3 r 2 = B, t r 2 T 3 , . (3.13) 

The (Ai) and (/!,) are linear combinations (using the g“)of the (f /t ) and the (m M ). We have 
normalized by setting r L = 1 and divided some of the equations by r 2 and r 3 . 

This gives 


/"■N 


i 4 2 

7i — ^1,72 = •—,73 
T 2 


7-3 


i?2 03 

fi = 0i, — —■, q — — 

7*2 7- 3 





_(A. t - 6Bi)r3 _ 

(Ai - 50i)r 3 + (/1 3 - 5B 3 ) . 


(3-14) 


where S = A 2 /B 2 . We can substitute this into (3.13) to obtain a quartic equation for r 3 
which will have at most four solutions. Then we solve (3.14) for r 2 and (3.11) and (3.12) for 
the ( 7 ,-) and (q). 

The basis matrices defined by (3.9) were choosen arbitrarily. The same analysis could be 
applied for other basis matrices. For some choices of basis there is only one possible 
solution. An example is 


'! 

0 

O' 


'0 

1 

O' 


'0 

0 

r 


'0 

0 

O' 

0 

0 

0 

> 

1 

0 

0 

, 

0 

0 

0 

} 

0 

0 

1 

.0 

0 

0. 


.0 

0 

0. 


.1 

0 

0. 


.0 

1 

0. 


This basis, however, requires that 

J = 0, J a' 1 0|(X)dX = 0. 


(3.15) 


(3.16) 


which seems an unnatural choice for the (IT). This illustrates the important point that the 
choice of basis functions is crucial. By assuming that we can express our illumination 
and surface reflectance in terms of a finite set of basis functions we are performing an 
implicit type of regularization (Poggio and Torre, 1984) by restricting the space of possible 
solutions. Only if the basis functions are choosen so that linear combinations of them are 
able to represent most of the naturally occuring illuminants and surface reflectances will 
methods of this type give good solutions. 
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This strategy of taking linear combinations of the ('/’T) will drastically simplify the equations. 
f"*'' It is also possible to solve equations (3.8) directly. We write the equation as 


= 7 i(T'nn + T'Vs + r? 3 r 3 ) + 7 8 (^ l r 1 + T£ 2 r 2 + It 3 r 3 ) + 73(^1 + ^ 2 r 2 + = 1, ..4 

( 3 - 17 ) 

We can use three of these four equations to solve for ( 7 ,) in terms of the (r<). This will leave 
us with two polynomial equations relating r 2 and r 3 (we have set n — 1). Similarly we solve 
similar equations for (/?,•) in terms of the (a,-) leaving another polynomial equation in r 2 and 
r 3 . These two equations in r 2 and r 3 can be combined to give solutions for t 2 and r 3 . Since 
these equations are polynomials there v/ill usually be more than one solution, and in some 
cases an infinite number. The previous method involving basis matrices is a lot simpler. If 
it is impossible to choose a basis like (3.9) it will still be possible to simplify the equations 
by taking linear combinations of the (7’T) to reduce the matrices to simple forms. 

Note that so far nothing has been said about the behaviour of the solution. It is even 
conceivable, though unlikely, that the solutions are sometimes not positive everywhere. 
They are guaranteed to be smooth provided the basis functions are. Instead of choosing 
three basis functions and solving for their coefficients it would be possible to have more 
basis functions and impose some a priori expectations on the solutions. Such an approach 
is discussed in Richards, Yuille and Poggio (1985). 

The Spatial Extension 

In this section we consider extensions to situations when the surface reflectance and the 
illumination are functions of space. Hurlbert and Poggio (1905) have shown that for many 
cases the spatial dependence of the illumination and surface reflectance can be factored 
out from the spectral dependence for colour patches. For these cases we can write 


S(x,\) = S(\)f{x) 


(4.1) 


and 


E[x, X) = E(\)g(x) 


(4.2) 


where f(x) and g(x) depend only on position. 

With these assumptions the spatial dependence of the illuminant and the surface reflectance, 
when not due to material (colour) changes, is factored out and the analysis can proceed 
as in the previous section. There is an important difference, however, the scaling ambiguity 
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(2.7) will now have a spatial dependence, In other words we will not be able to distinguish 
the spatial dependences f(x) and g[x) of the surface reflectance and illumination. For 
many cases the illumination will have little spatial variation (for example, if the light sources 
are a long way away) and we can set g(x) = 1 and solve for f(x). Otherwise we can 
consider the boundaries between two materials and assume that the spatial dependence 
of the illumination is negligible across the boundary. We can then solve for the colour on 
either side of the boundaries (and fill in if necessary). A third method is to note that if there 
is a single light source the spatial dependence g(x) will have a definite form , typically it will 
fall off as 1/r 2 (where r is the distance to the light source), and so f(x) can be determined. 

In the previous sections we were considering a Mondrian style world consisting of flat 
patches of uniform colour. Our method would sample two patches and obtain two equations 
(2.5) and (2.6). For the more general case we would sample neighbouring points. If the 
two points lie on the same material the l„ will be proportional to the and the equations 
will be underdetermined. If the points occur at a material boundary, one lying in each 
material, then we have enough equations to solve for the colours. Since we are considering 
neighbouring points it is reasonable to assume that the illuminant is the same (including the 
scaling factor) for each point. 

To summarize: we sample neighbouring points in the image and compare their values 
and If these are proportional we conclude that the points have the same colour (though 
not necessarily the same brightness) and hence lie on the same material. If they are not 
then we are on a material boundary and can solve for the colours on either sides of the 
boundary. The colour cannot change inside the boundary, or else we would detect it, and so 
the colour on the boundaries determines the colour in the interior. However the brightness, 
or strength, of the colour is undetermined. If the illuminant is constant in space we can 
solve this by a single global normalization (such as setting n = 1, as in (3.5)). We suggest 
making this assumption even if the illuminant varies spatially, unless the variation is very 
large or there is strong evidence to the contrary from other souces. Thus spatial variation 
of the illuminant will be interpreted as due to changes in the reflectance function due to 
orientation and other changes. 

Other ways of generalizing to allow for spatial dependence would include using the standard 
assumption that the spatial variation of the illuminant is smaller than the spatial variation of 
the surface reflectance (Land (1983), Horn (1974), Blake (1984)). New methods for doing 
this are discussed in Hurlbert and Poggio (1985) and extensions will be discussed in a 
forthcoming paper. 
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